require(cjoint)
require(lfe)
require(aod)
require(RColorBrewer)
Set1 <- brewer.pal(9, "Set1")
windowsFonts(Helvetica = windowsFont("Helvetica"))

attribute.label <- c("Parents", "Identity", "Activity", 
                     "Policy", "Habit", "Intention")
level.label <- list(c("Independent", "Same party", "Opponent party"), 
                    c("Does not matter", "Personal insult"), 
                    c("Never", "Rarely", "Sometimes", "Often", "Very often"), 
                    c("Does not approve", "Partially approves", "Completely approves"), 
                    c("Always abstained", "Sometimes abstained", 
                      "Sometimes deviated", "Always loyal"), 
                    c("Abstain", "Same party", "Opponent party"))

choice.formula <- formula(choice ~ a.parents + b.identity + c.activity + d.policy + e.habit + f.intention)
rating.formula <- formula(rating ~ a.parents + b.identity + c.activity + d.policy + e.habit + f.intention)

#### load dataset ####
load("entire_data.Rdata")
load("clean_data.Rdata")

# exclude Independents
entire.partisan.data <- subset(entire.data, PID != 3)
partisan.data <- subset(clean.data, PID != 3)

#### original results for comparison ####
AMCE.choice.total.result <- summary(cjoint::amce(choice.formula, data = partisan.data, 
                                                 respondent.id = "respondent.id"))
AMCE.rating.total.result <- summary(cjoint::amce(rating.formula, data = partisan.data, 
                                                 respondent.id = "respondent.id"))

#### analysis using the data of the first round ####
AMCE.choice.first.result <- summary(cjoint::amce(choice.formula, 
                                                 data = partisan.data, 
                                                 subset = partisan.data$candidate.id < 3, 
                                                 respondent.id = "respondent.id"))
AMCE.rating.first.result <- summary(cjoint::amce(rating.formula, 
                                                 data = partisan.data, 
                                                 subset = partisan.data$candidate.id < 3, 
                                                 respondent.id = "respondent.id"))
# number of responses
AMCE.choice.first.result$samplesize_estimates

cairo_pdf("Figure_A4.pdf", width = 6.2, height = 3.8, pointsize = 7, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0), xpd = TRUE)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 5, 2, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
  points(0, plot.location + 0.15, pch = 19, col = Set1[1])
  points(0, plot.location - 0.15, pch = 15, col = Set1[2])
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
    points(AMCE.choice.total.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[1])
    segments(AMCE.choice.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.choice.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, col = Set1[1])
    points(AMCE.choice.first.result$amce$Estimate[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[2])
    segments(AMCE.choice.first.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.first.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, 
             AMCE.choice.first.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.first.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend(0.4, 26, legend = c("All rounds", "First task"), 
       pch = c(19, 15), col = Set1[1:2], bty = "n", xjust = 1, yjust = 0)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 5, 2, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
  points(0, plot.location + 0.15, pch = 19, col = Set1[1])
  points(0, plot.location - 0.15, pch = 15, col = Set1[2])
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
    points(AMCE.rating.total.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[1])
    segments(AMCE.rating.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.rating.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, col = Set1[1])
    points(AMCE.rating.first.result$amce$Estimate[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[2])
    segments(AMCE.rating.first.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.first.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, 
             AMCE.rating.first.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.first.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend(0.8, 26, legend = c("All rounds", "First task"), 
       pch = c(19, 15), col = Set1[1:2], bty = "n", xjust = 1, yjust = 0)
dev.off()

#### analysis including the fixed effects of the names of hypothetical persons ####
D.choice.result.wo.name <- felm(choice ~ a.parents + b.identity + c.activity + d.policy + 
                                  e.habit + f.intention | 0 | 0 | respondent.id, 
                                data = partisan.data, subset = partisan.data$condition == 0)
D.choice.result.w.name <- felm(choice ~ a.parents + b.identity + c.activity + d.policy + 
                                 e.habit + f.intention + candidate.name | 0 | 0 | respondent.id, 
                               data = partisan.data, subset = partisan.data$condition == 0)
D.rating.result.wo.name <- felm(rating ~ a.parents + b.identity + c.activity + d.policy + 
                                  e.habit + f.intention | 0 | 0 | respondent.id, 
                                data = partisan.data, subset = partisan.data$condition == 0)
D.rating.result.w.name <- felm(rating ~ a.parents + b.identity + c.activity + d.policy + 
                                 e.habit + f.intention + candidate.name | 0 | 0 | respondent.id, 
                               data = partisan.data, subset = partisan.data$condition == 0)
R.choice.result.wo.name <- felm(choice ~ a.parents + b.identity + c.activity + d.policy + 
                                  e.habit + f.intention | 0 | 0 | respondent.id, 
                                data = partisan.data, subset = partisan.data$condition == 1)
R.choice.result.w.name <- felm(choice ~ a.parents + b.identity + c.activity + d.policy + 
                                 e.habit + f.intention + candidate.name | 0 | 0 | respondent.id, 
                               data = partisan.data, subset = partisan.data$condition == 1)
R.rating.result.wo.name <- felm(rating ~ a.parents + b.identity + c.activity + d.policy + 
                                  e.habit + f.intention | 0 | 0 | respondent.id, 
                                data = partisan.data, subset = partisan.data$condition == 1)
R.rating.result.w.name <- felm(rating ~ a.parents + b.identity + c.activity + d.policy + 
                                 e.habit + f.intention + candidate.name | 0 | 0 | respondent.id, 
                               data = partisan.data, subset = partisan.data$condition == 1)

cairo_pdf("Figure_A5.pdf", width = 6.4, height = 2, pointsize = 6, family = "Helvetica")
layout(matrix(1:4, 1, 4))
par(mar = c(2, 0, 3, 0), xpd = TRUE)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice\n(D Condition)", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
  points(0, plot.location + 0.15, pch = 19, col = Set1[1])
  points(0, plot.location - 0.15, pch = 15, col = Set1[2])
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
    points(D.choice.result.wo.name$coefficients[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[1])
    segments(D.choice.result.wo.name$coefficients[coefficients.location] + 
               qnorm(0.025) * D.choice.result.wo.name$cse[coefficients.location], 
             plot.location + 0.15, 
             D.choice.result.wo.name$coefficients[coefficients.location] + 
               qnorm(0.975) * D.choice.result.wo.name$cse[coefficients.location], 
             plot.location + 0.15, col = Set1[1])
    points(D.choice.result.w.name$coefficients[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[2])
    segments(D.choice.result.w.name$coefficients[coefficients.location] + 
               qnorm(0.025) * D.choice.result.w.name$cse[coefficients.location], 
             plot.location - 0.15, 
             D.choice.result.w.name$coefficients[coefficients.location] + 
               qnorm(0.975) * D.choice.result.w.name$cse[coefficients.location], 
             plot.location - 0.15, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend(0.4, 26, legend = c("w/o name FE", "w/ name FE"), 
       pch = c(19, 15), col = Set1[1:2], bty = "n", xjust = 1, yjust = 0)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating\n(D Condition)", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
  points(0, plot.location + 0.15, pch = 19, col = Set1[1])
  points(0, plot.location - 0.15, pch = 15, col = Set1[2])
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
    points(D.rating.result.wo.name$coefficients[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[1])
    segments(D.rating.result.wo.name$coefficients[coefficients.location] + 
               qnorm(0.025) * D.rating.result.wo.name$cse[coefficients.location], 
             plot.location + 0.15, 
             D.rating.result.wo.name$coefficients[coefficients.location] + 
               qnorm(0.975) * D.rating.result.wo.name$cse[coefficients.location], 
             plot.location + 0.15, col = Set1[1])
    points(D.rating.result.w.name$coefficients[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[2])
    segments(D.rating.result.w.name$coefficients[coefficients.location] + 
               qnorm(0.025) * D.rating.result.w.name$cse[coefficients.location], 
             plot.location - 0.15, 
             D.rating.result.w.name$coefficients[coefficients.location] + 
               qnorm(0.975) * D.rating.result.w.name$cse[coefficients.location], 
             plot.location - 0.15, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend(0.8, 26, legend = c("w/o name FE", "w/ name FE"), 
       pch = c(19, 15), col = Set1[1:2], bty = "n", xjust = 1, yjust = 0)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice\n(R Condition)", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
  points(0, plot.location + 0.15, pch = 19, col = Set1[1])
  points(0, plot.location - 0.15, pch = 15, col = Set1[2])
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
    points(R.choice.result.wo.name$coefficients[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[1])
    segments(R.choice.result.wo.name$coefficients[coefficients.location] + 
               qnorm(0.025) * R.choice.result.wo.name$cse[coefficients.location], 
             plot.location + 0.15, 
             R.choice.result.wo.name$coefficients[coefficients.location] + 
               qnorm(0.975) * R.choice.result.wo.name$cse[coefficients.location], 
             plot.location + 0.15, col = Set1[1])
    points(R.choice.result.w.name$coefficients[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[2])
    segments(R.choice.result.w.name$coefficients[coefficients.location] + 
               qnorm(0.025) * R.choice.result.w.name$cse[coefficients.location], 
             plot.location - 0.15, 
             R.choice.result.w.name$coefficients[coefficients.location] + 
               qnorm(0.975) * R.choice.result.w.name$cse[coefficients.location], 
             plot.location - 0.15, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend(0.4, 26, legend = c("w/o name FE", "w/ name FE"), 
       pch = c(19, 15), col = Set1[1:2], bty = "n", xjust = 1, yjust = 0)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating\n(R Condition)", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
  points(0, plot.location + 0.15, pch = 19, col = Set1[1])
  points(0, plot.location - 0.15, pch = 15, col = Set1[2])
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
    points(R.rating.result.wo.name$coefficients[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[1])
    segments(R.rating.result.wo.name$coefficients[coefficients.location] + 
               qnorm(0.025) * R.rating.result.wo.name$cse[coefficients.location], 
             plot.location + 0.15, 
             R.rating.result.wo.name$coefficients[coefficients.location] + 
               qnorm(0.975) * R.rating.result.wo.name$cse[coefficients.location], 
             plot.location + 0.15, col = Set1[1])
    points(R.rating.result.w.name$coefficients[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[2])
    segments(R.rating.result.w.name$coefficients[coefficients.location] + 
               qnorm(0.025) * R.rating.result.w.name$cse[coefficients.location], 
             plot.location - 0.15, 
             R.rating.result.w.name$coefficients[coefficients.location] + 
               qnorm(0.975) * R.rating.result.w.name$cse[coefficients.location], 
             plot.location - 0.15, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend(0.8, 26, legend = c("w/o name FE", "w/ name FE"), 
       pch = c(19, 15), col = Set1[1:2], bty = "n", xjust = 1, yjust = 0)
dev.off()

#### analysis not discarding satisficing responses ####
AMCE.choice.all.result <- summary(cjoint::amce(choice.formula, 
                                               data = entire.partisan.data, 
                                               respondent.id = "respondent.id"))
AMCE.rating.all.result <- summary(cjoint::amce(rating.formula, 
                                               data = entire.partisan.data, 
                                               respondent.id = "respondent.id"))

# numbers of responses and respondents
nrow(entire.partisan.data)
length(unique(entire.partisan.data$respondent.id))

cairo_pdf("Figure_A6.pdf", width = 6.2, height = 3.8, pointsize = 7, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0), xpd = TRUE)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
  points(0, plot.location + 0.15, pch = 19, col = Set1[1])
  points(0, plot.location - 0.15, pch = 15, col = Set1[2])
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
    points(AMCE.choice.total.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[1])
    segments(AMCE.choice.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.choice.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, col = Set1[1])
    points(AMCE.choice.all.result$amce$Estimate[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[2])
    segments(AMCE.choice.all.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.all.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, 
             AMCE.choice.all.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.all.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend(0.4, 26, legend = c("w/o satisficing", "w/ satisficing"), 
       pch = c(19, 15), col = Set1[1:2], bty = "n", xjust = 1, yjust = 0)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
  points(0, plot.location + 0.15, pch = 19, col = Set1[1])
  points(0, plot.location - 0.15, pch = 15, col = Set1[2])
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
    points(AMCE.rating.total.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[1])
    segments(AMCE.rating.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.rating.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, col = Set1[1])
    points(AMCE.rating.all.result$amce$Estimate[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[2])
    segments(AMCE.rating.all.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.all.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, 
             AMCE.rating.all.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.all.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend(0.8, 26, legend = c("w/o satisficing", "w/ satisficing"), 
       pch = c(19, 15), col = Set1[1:2], bty = "n", xjust = 1, yjust = 0)
dev.off()

#### analysis discarding respondents who failed to pass the IMC ####
# discard respondents who failed to pass the IMC
screened.partisan.data <- subset(entire.partisan.data, 
                                 directed == TRUE & IMC == TRUE & time == TRUE)

# number of responses
nrow(screened.partisan.data)
round(nrow(screened.partisan.data) / nrow(entire.partisan.data), 3) * 100
# number of respondents
length(unique(screened.partisan.data$respondent.id))
round(length(unique(screened.partisan.data$respondent.id)) / 
        length(unique(entire.partisan.data$respondent.id)), 3)

AMCE.choice.perfect.result <- summary(cjoint::amce(choice.formula, 
                                                   data = screened.partisan.data, 
                                                   respondent.id = "respondent.id"))
AMCE.rating.perfect.result <- summary(cjoint::amce(rating.formula, 
                                                   data = screened.partisan.data, 
                                                   respondent.id = "respondent.id"))

cairo_pdf("Figure_A7.pdf", width = 6.2, height = 3.8, pointsize = 7, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0), xpd = TRUE)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
  points(0, plot.location + 0.15, pch = 19, col = Set1[1])
  points(0, plot.location - 0.15, pch = 15, col = Set1[2])
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
    points(AMCE.choice.total.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[1])
    segments(AMCE.choice.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.choice.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, col = Set1[1])
    points(AMCE.choice.perfect.result$amce$Estimate[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[2])
    segments(AMCE.choice.perfect.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.perfect.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, 
             AMCE.choice.perfect.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.perfect.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend(0.4, 26, legend = c("Not using IMC", "Using IMC"), 
       pch = c(19, 15), col = Set1[1:2], bty = "n", xjust = 1, yjust = 0)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
  points(0, plot.location + 0.15, pch = 19, col = Set1[1])
  points(0, plot.location - 0.15, pch = 15, col = Set1[2])
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9])
    points(AMCE.rating.total.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[1])
    segments(AMCE.rating.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.rating.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, col = Set1[1])
    points(AMCE.rating.perfect.result$amce$Estimate[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[2])
    segments(AMCE.rating.perfect.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.perfect.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, 
             AMCE.rating.perfect.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.perfect.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend(0.8, 26, legend = c("Not using IMC", "Using IMC"), 
       pch = c(19, 15), col = Set1[1:2], bty = "n", xjust = 1, yjust = 0)
dev.off()

#### analysis of interaction effects between multiple attributes ####
# function for joint-significance tests
interaction.test <- function(variable.1, variable.2, data = data, choice = TRUE) {
  nfactor.1 <- nlevels(eval(parse(text = paste0("partisan.data$", variable.1))))
  nfactor.2 <- nlevels(eval(parse(text = paste0("partisan.data$", variable.2))))
  if (choice == TRUE) {
    felm.formula <- as.formula(paste0("choice ~ a.parents + b.identity + c.activity + 
                                      d.policy + e.habit + f.intention + ", 
                                      variable.1, " : ", variable.2, 
                                      " | 0 | 0 | respondent.id"))
  } else {
    felm.formula <- as.formula(paste0("rating ~ a.parents + b.identity + c.activity + 
                                      d.policy + e.habit + f.intention + ", 
                                      variable.1, " : ", variable.2, 
                                      " | 0 | 0 | respondent.id"))
  }
  felm.result <- felm(felm.formula, data = data)
  Wald.test <- wald.test(Sigma = felm.result$clustervcv, 
                         b = felm.result$coefficients, 
                         Terms = (length(felm.result$coefficients) - 
                                    (nfactor.1 - 1) * (nfactor.2 - 1) + 1):
                           length(felm.result$coefficients))
  Wald.test$result$chi2[3]
}

variable.label <- c("a.parents", "b.identity", "c.activity", 
                    "d.policy", "e.habit", "f.intention")
interaction.test.p.choice <- interaction.test.p.rating <- matrix(NA, 6, 6)
for (i in 1:6) {
  for (j in 1:6) {
    if (i == j) break
    interaction.test.p.choice[i, j] <- interaction.test(variable.label[i], 
                                                        variable.label[j], 
                                                        data = partisan.data)
    interaction.test.p.rating[i, j] <- interaction.test(variable.label[i], 
                                                        variable.label[j], 
                                                        data = partisan.data, 
                                                        choice = FALSE)
  }
}

# Table A.1
round(interaction.test.p.choice, 3)
round(interaction.test.p.rating, 3)

#### Analysis of self-identified Independents ####
AMCE.choice.partisan.result <- summary(cjoint::amce(choice.formula, 
                                                    data = clean.data, 
                                                    subset = clean.data$PID != 3, 
                                                    respondent.id = "respondent.id"))
AMCE.rating.partisan.result <- summary(cjoint::amce(rating.formula, 
                                                    data = clean.data, 
                                                    subset = clean.data$PID != 3, 
                                                    respondent.id = "respondent.id"))
AMCE.choice.Ind.result <- summary(cjoint::amce(choice.formula, 
                                               data = clean.data, 
                                               subset = clean.data$PID == 3, 
                                               respondent.id = "respondent.id"))
AMCE.rating.Ind.result <- summary(cjoint::amce(rating.formula, 
                                               data = clean.data, 
                                               subset = clean.data$PID == 3, 
                                               respondent.id = "respondent.id"))

# number of responses
AMCE.choice.total.result$samplesize_estimates
AMCE.choice.Ind.result$samplesize_estimates
# number of respondents
AMCE.choice.total.result$respondents
AMCE.choice.Ind.result$respondents

cairo_pdf("Figure_A8.pdf", width = 6.2, height = 3.8, pointsize = 7, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0), xpd = TRUE)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location + 0.15, pch = 19)
  points(0, plot.location - 0.15, pch = 15)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(AMCE.choice.total.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19)
    segments(AMCE.choice.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.choice.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15)
    points(AMCE.choice.Ind.result$amce$Estimate[coefficients.location], plot.location - 0.15, 
           pch = 15)
    segments(AMCE.choice.Ind.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.Ind.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, 
             AMCE.choice.Ind.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.Ind.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend(0.4, 26, legend = c("Democrat and Republican", "Independent"), 
       pch = c(19, 15), col = c(1, 1), bty = "n", xjust = 1, yjust = 0, cex = 0.8)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location + 0.15, pch = 19)
  points(0, plot.location - 0.15, pch = 15)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(AMCE.rating.total.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19)
    segments(AMCE.rating.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.rating.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15)
    points(AMCE.rating.Ind.result$amce$Estimate[coefficients.location], plot.location - 0.15, 
           pch = 15)
    segments(AMCE.rating.Ind.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.Ind.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, 
             AMCE.rating.Ind.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.Ind.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend(0.8, 26, legend = c("Democrat and Republican", "Independent"), 
       pch = c(19, 15), col = c(1, 1), bty = "n", xjust = 1, yjust = 0, cex = 0.8)
dev.off()

#### analysis with calibration weights ####
clean.data.svy.acs <- svydesign(ids = ~ respondent.id, 
                                data = clean.data[! is.na(clean.data$weight.acs), ],
                                weights = ~ weight.acs)
AMCE.choice.acs <- summary(svyglm(choice.formula, 
                                  design = clean.data.svy.acs))
AMCE.rating.acs <- summary(svyglm(rating.formula, 
                                  design = clean.data.svy.acs))

clean.data.svy.cps <- svydesign(ids = ~ respondent.id, 
                                data = clean.data[! is.na(clean.data$weight.cps), ],
                                weights = ~ weight.cps)
AMCE.choice.cps <- summary(svyglm(choice.formula, 
                                  design = clean.data.svy.cps))
AMCE.rating.cps <- summary(svyglm(rating.formula, 
                                  design = clean.data.svy.cps))

# results using the entire sample including Independents for comparison
AMCE.choice.entire.result <- summary(cjoint::amce(choice.formula, data = clean.data, 
                                                  respondent.id = "respondent.id"))
AMCE.rating.entire.result <- summary(cjoint::amce(rating.formula, data = clean.data, 
                                                  respondent.id = "respondent.id"))

cairo_pdf("Figure_A9.pdf", width = 6.2, height = 3.8, pointsize = 7, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0))
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location + 0.15, pch = 19, cex = 0.7)
  points(0, plot.location, pch = 15, col = Set1[1], cex = 0.7)
  points(0, plot.location - 0.15, pch = 17, col = Set1[2], cex = 0.7)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(AMCE.choice.entire.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, cex = 0.7)
    segments(AMCE.choice.entire.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.entire.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.choice.entire.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.entire.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, lwd = 0.7)
    points(AMCE.choice.acs$coefficients[coefficients.location + 1, 1], plot.location, 
           pch = 15, col = Set1[1], cex = 0.7)
    segments(AMCE.choice.acs$coefficients[coefficients.location + 1, 1] + 
               qnorm(0.025) * AMCE.choice.acs$coefficients[coefficients.location + 1, 2], 
             plot.location, 
             AMCE.choice.acs$coefficients[coefficients.location + 1, 1] + 
               qnorm(0.975) * AMCE.choice.acs$coefficients[coefficients.location + 1, 2], 
             plot.location, lwd = 0.7, col = Set1[1])
    points(AMCE.choice.cps$coefficients[coefficients.location + 1, 1], plot.location - 0.15, 
           pch = 17, col = Set1[2], cex = 0.7)
    segments(AMCE.choice.cps$coefficients[coefficients.location + 1, 1] + 
               qnorm(0.025) * AMCE.choice.cps$coefficients[coefficients.location + 1, 2], 
             plot.location - 0.15, 
             AMCE.choice.cps$coefficients[coefficients.location + 1, 1] + 
               qnorm(0.975) * AMCE.choice.cps$coefficients[coefficients.location + 1, 2], 
             plot.location - 0.15, lwd = 0.7, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend("topright", legend = c("w/o weights", "ACS weights", "CPS weights"), 
       pch = c(19, 15, 17), col = c(1, Set1[1], Set1[2]), bty = "n", xjust = 1, yjust = 0, cex = 0.6)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location + 0.15, pch = 19, cex = 0.7)
  points(0, plot.location, pch = 15, col = Set1[1], cex = 0.7)
  points(0, plot.location - 0.15, pch = 17, col = Set1[2], cex = 0.7)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(AMCE.rating.entire.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, cex = 0.7)
    segments(AMCE.rating.entire.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.entire.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.rating.entire.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.entire.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, lwd = 0.7)
    points(AMCE.rating.acs$coefficients[coefficients.location + 1, 1], plot.location, 
           pch = 15, col = Set1[1], cex = 0.7)
    segments(AMCE.rating.acs$coefficients[coefficients.location + 1, 1] + 
               qnorm(0.025) * AMCE.rating.acs$coefficients[coefficients.location + 1, 2], 
             plot.location, 
             AMCE.rating.acs$coefficients[coefficients.location + 1, 1] + 
               qnorm(0.975) * AMCE.rating.acs$coefficients[coefficients.location + 1, 2], 
             plot.location, lwd = 0.7, col = Set1[1])
    points(AMCE.rating.cps$coefficients[coefficients.location + 1, 1], plot.location - 0.15, 
           pch = 17, col = Set1[2], cex = 0.7)
    segments(AMCE.rating.cps$coefficients[coefficients.location + 1, 1] + 
               qnorm(0.025) * AMCE.rating.cps$coefficients[coefficients.location + 1, 2], 
             plot.location - 0.15, 
             AMCE.rating.cps$coefficients[coefficients.location + 1, 1] + 
               qnorm(0.975) * AMCE.rating.cps$coefficients[coefficients.location + 1, 2], 
             plot.location - 0.15, lwd = 0.7, col = Set1[2])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend("topright", legend = c("w/o weights", "ACS weights", "CPS weights"), 
       pch = c(19, 15, 17), col = c(1, Set1[1], Set1[2]), bty = "n", xjust = 1, yjust = 0, cex = 0.6)
dev.off()